Load packages required for analysis
Load the fitted model and prior draws.
if (class(model) %in% "character") {
model <- ModelTBBCGEngland::read_libbi(file.path(model_dir, "libbi", "posterior"))
}
priors <- readRDS(file.path(model_dir, "data", "prior-params.rds"))
model
## Wrapper around LibBi
## ======================
## Model: Baseline
## Run time: 406.828 seconds
## Number of samples: 1000
## State trajectories recorded: E H L P S T_E T_P YearlyAgeCases YearlyCases YearlyEAgeCases YearlyECases YearlyEPulCases YearlyEPulDeaths YearlyNonUKborn YearlyPAgeCases YearlyPCases YearlyPulCases YearlyPulDeaths
## Observation trajectories recorded: YearlyAgeInc YearlyHistPInc YearlyInc
## Parameters recorded: alpha_t_decay alpha_t_init c_eff c_hist chi_init delta epsilon_h_0_4 epsilon_h_15_89 epsilon_h_5_14 epsilon_l_0_4 epsilon_l_15_89 epsilon_l_5_14 HistMeasError kappa_0_4 kappa_15_89 kappa_5_14 M mu_e_0_14 mu_e_15_59 mu_e_60_89 mu_p_0_14 mu_p_15_59 mu_p_60_89 nu_e_0_14 nu_e_15_89 nu_p_0_14 nu_p_15_89 phi_0_14 phi_15_59 phi_60_89 rho_0_14 rho_15_59 rho_60_89 Upsilon_0_14 Upsilon_15_59 Upsilon_60_89 zeta_e_0_14 zeta_e_15_59 zeta_e_60_89 zeta_p_0_14 zeta_p_15_59 zeta_p_60_89
traces <- coda::mcmc(rbi::get_traces(model))
coda::rejectionRate(traces)
## M c_eff c_hist delta
## 0.995996 0.995996 0.995996 0.995996
## epsilon_h_0_4 epsilon_h_5_14 epsilon_h_15_89 kappa_0_4
## 0.995996 0.995996 0.995996 0.995996
## kappa_5_14 kappa_15_89 epsilon_l_0_4 epsilon_l_5_14
## 0.995996 0.995996 0.995996 0.995996
## epsilon_l_15_89 phi_0_14 phi_15_59 phi_60_89
## 0.995996 0.995996 0.995996 0.995996
## Upsilon_0_14 Upsilon_15_59 Upsilon_60_89 rho_0_14
## 0.995996 0.995996 0.995996 0.995996
## rho_15_59 rho_60_89 nu_p_0_14 nu_p_15_89
## 0.995996 0.995996 0.995996 0.995996
## nu_e_0_14 nu_e_15_89 zeta_p_0_14 zeta_p_15_59
## 0.995996 0.995996 0.995996 0.995996
## zeta_p_60_89 zeta_e_0_14 zeta_e_15_59 zeta_e_60_89
## 0.995996 0.995996 0.995996 0.995996
## mu_p_0_14 mu_p_15_59 mu_p_60_89 mu_e_0_14
## 0.995996 0.995996 0.995996 0.995996
## mu_e_15_59 mu_e_60_89 chi_init alpha_t_init
## 0.995996 0.995996 0.995996 0.000000
## alpha_t_decay HistMeasError
## 0.995996 0.995996
coda::autocorr.plot(traces)
- Evaluate Posterior Traces
plot_param(model, scales = "free", plot_type = "trace", burn_in = 0) + theme(legend.position = "none")
plot_param(model, prior_params = priors, scales = "free")
param_sum <- plot_param(model, prior_params = priors, plot_data = FALSE) %>%
group_by(Distribution, parameter, length) %>%
summarise(median = median(value),
lll = quantile(value, 0.025),
hhh = quantile(value, 0.975)) %>%
mutate(value = pretty_ci(median, lll, hhh, sep = ", ")) %>%
group_by(Distribution, parameter) %>%
mutate(Parameter = case_when(max(length) > 1 ~ paste(parameter, 1:n(), sep = "-"),
TRUE ~ as.character(parameter))
) %>%
ungroup %>%
select(Distribution, Parameter, value) %>%
spread(key = "Distribution", value = "value") %>%
select(Parameter, Prior, Posterior)
saveRDS(param_sum, file.path(model_dir, "data", "sum_prior_posterior.rds"))
knitr::kable(param_sum)
| Parameter | Prior | Posterior |
|---|---|---|
| alpha_t_decay | -0.13 (-0.24, -0.03) | -0.19 (-0.26, -0.05) |
| alpha_t_init | 0.84 (0.76, 0.90) | 0.50 (0.11, 0.89) |
| c_eff | 2.46 (0.14, 4.88) | 3.84 (1.43, 4.18) |
| c_hist | 12.58 (10.12, 14.89) | 12.00 (10.81, 13.45) |
| chi_init | 0.19 (0.08, 0.29) | 0.23 (0.18, 0.30) |
| delta | 0.78 (0.70, 0.85) | 0.78 (0.77, 0.84) |
| epsilon_h_0_4 | 1.53 (0.10, 4.88) | 1.11 (0.69, 3.72) |
| epsilon_h_15_89 | 25.72 (1.45, 77.21) | 9.04 (4.43, 70.81) |
| epsilon_h_5_14 | 3.72 (0.16, 11.99) | 5.85 (4.21, 7.45) |
| epsilon_l_0_4 | 594.41 (29.25, 1714.03) | 325.02 (183.07, 1257.01) |
| epsilon_l_15_89 | 1066.55 (38.46, 3368.79) | 44.09 (31.10, 51.12) |
| epsilon_l_5_14 | 537.07 (31.69, 1541.71) | 837.75 (28.05, 1249.42) |
| HistMeasError | 1.00 (0.60, 1.37) | 1.16 (1.03, 1.38) |
| kappa_0_4 | 0.88 (0.05, 2.78) | 1.01 (0.06, 1.21) |
| kappa_15_89 | 1.11 (0.04, 3.47) | 1.26 (0.62, 1.87) |
| kappa_5_14 | 0.97 (0.05, 3.10) | 0.61 (0.18, 0.73) |
| M | 0.25 (0.01, 0.49) | 0.21 (0.01, 0.40) |
| mu_e_0_14 | 275.27 (211.92, 336.54) | 256.76 (249.02, 301.04) |
| mu_e_15_59 | 194.13 (69.41, 321.66) | 160.81 (127.82, 282.82) |
| mu_e_60_89 | 36.50 (1.43, 106.56) | 35.09 (0.80, 42.95) |
| mu_p_0_14 | 240.25 (154.15, 330.24) | 248.58 (195.58, 270.37) |
| mu_p_15_59 | 88.24 (5.27, 259.49) | 41.89 (14.04, 68.54) |
| mu_p_60_89 | 46.99 (2.35, 148.94) | 75.42 (41.67, 127.12) |
| nu_e_0_14 | 0.17 (0.00, 1.30) | 0.15 (0.04, 0.26) |
| nu_e_15_89 | 0.33 (0.01, 1.67) | 0.71 (0.08, 2.72) |
| nu_p_0_14 | 0.13 (0.00, 0.71) | 0.07 (0.01, 0.36) |
| nu_p_15_89 | 0.23 (0.01, 1.17) | 0.16 (0.06, 0.52) |
| phi_0_14 | 0.58 (0.30, 1.06) | 0.47 (0.40, 0.55) |
| phi_15_59 | 0.62 (0.30, 1.15) | 0.38 (0.20, 0.78) |
| phi_60_89 | 0.60 (0.29, 1.12) | 0.68 (0.31, 0.97) |
| rho_0_14 | 0.30 (0.26, 0.34) | 0.29 (0.27, 0.33) |
| rho_15_59 | 0.65 (0.64, 0.66) | 0.65 (0.64, 0.65) |
| rho_60_89 | 0.54 (0.52, 0.55) | 0.53 (0.53, 0.54) |
| Upsilon_0_14 | 0.63 (0.61, 0.65) | 0.64 (0.62, 0.64) |
| Upsilon_15_59 | 0.71 (0.70, 0.71) | 0.71 (0.70, 0.71) |
| Upsilon_60_89 | 0.75 (0.74, 0.76) | 0.75 (0.75, 0.76) |
| zeta_e_0_14 | 123.74 (58.77, 187.11) | 174.99 (106.52, 193.46) |
| zeta_e_15_59 | 62.71 (3.84, 175.77) | 88.10 (5.55, 156.23) |
| zeta_e_60_89 | 132.08 (55.59, 209.75) | 90.31 (84.13, 191.14) |
| zeta_p_0_14 | 97.26 (17.74, 184.18) | 89.02 (51.59, 107.26) |
| zeta_p_15_59 | 81.34 (3.53, 252.31) | 133.28 (27.51, 146.85) |
| zeta_p_60_89 | 122.28 (16.92, 246.71) | 155.70 (10.61, 231.57) |
plot_state(model, summarise = TRUE)
plot_state(model, summarise = TRUE, start_time = 59)
plot_state(model, states = c("YearlyHistPInc", "YearlyInc"), summarise = TRUE, start_time = 49) -> p2
p2
pred_obs <- plot_state(model, states = c("YearlyHistPInc", "YearlyInc", "YearlyAgeInc"), summarise = FALSE, start_time = 59, plot_data = FALSE)
obs <- ModelTBBCGEngland::setup_model_obs() %>%
bind_rows(.id = "state") %>%
mutate(time = time + 1931)
pred_obs <- pred_obs %>%
dplyr::filter(Average == "median") %>%
mutate(pred_value = pretty_ci(Count, lll, hhh, digits = 0)) %>%
left_join(obs, by = c("state" = "state", "time" = "time", "age" = "age")) %>%
select(Observation = state, Age = age, Year = time, `Observed Incidence` = value, `Predicted Incidence` = pred_value) %>%
mutate(Year = Year)
sum_obs_table <- pred_obs %>%
dplyr::filter(Observation %in% c("YearlyInc", "YearlyHistPInc")) %>%
drop_na(`Observed Incidence`) %>%
select(-Age) %>%
mutate(Observation = case_when(Observation == "YearlyInc" ~ "All TB cases",
Observation == "YearlyHistPInc" ~ "Pulmonary TB cases"))
saveRDS(sum_obs_table, file.path(model_dir, "data", "sum_obs_table.rds"))
age_cases_table <- pred_obs %>%
dplyr::filter(Observation %in% c("YearlyAgeInc")) %>%
drop_na(`Observed Incidence`) %>%
select(-Observation) %>%
group_by(Year) %>%
mutate(Age = c(paste0(seq(0, 65, 5), "-", seq(4, 69, 5)), "70-89")) %>%
ungroup
saveRDS(age_cases_table, file.path(model_dir, "data", "age_cases_table.rds"))
knitr::kable(sum_obs_table)
| Observation | Year | Observed Incidence | Predicted Incidence |
|---|---|---|---|
| Pulmonary TB cases | 1990 | 3618 | 4181 (2798 to 4349) |
| Pulmonary TB cases | 1991 | 3596 | 4100 (2745 to 4400) |
| Pulmonary TB cases | 1992 | 3816 | 3905 (2803 to 4318) |
| Pulmonary TB cases | 1993 | 3961 | 3790 (2604 to 4275) |
| Pulmonary TB cases | 1994 | 3694 | 3731 (2558 to 4103) |
| Pulmonary TB cases | 1995 | 3711 | 3709 (2614 to 4067) |
| Pulmonary TB cases | 1996 | 3690 | 3437 (2451 to 4033) |
| Pulmonary TB cases | 1997 | 3947 | 3503 (2431 to 3988) |
| Pulmonary TB cases | 1998 | 3926 | 3304 (2434 to 3910) |
| Pulmonary TB cases | 1999 | 3827 | 3286 (2384 to 3807) |
| All TB cases | 2000 | 1803 | 2417 (1618 to 3261) |
| All TB cases | 2001 | 1866 | 2319 (1521 to 3107) |
| All TB cases | 2002 | 1833 | 2105 (1521 to 3102) |
| All TB cases | 2003 | 1685 | 2124 (1371 to 3090) |
| All TB cases | 2004 | 1776 | 1949 (1290 to 3058) |
knitr::kable(age_cases_table)
| Age | Year | Observed Incidence | Predicted Incidence |
|---|---|---|---|
| 0-4 | 2000 | 75 | 1 (0 to 14) |
| 5-9 | 2000 | 59 | 0 (0 to 12) |
| 10-14 | 2000 | 75 | 1 (0 to 21) |
| 15-19 | 2000 | 112 | 4 (2 to 35) |
| 20-24 | 2000 | 133 | 6 (2 to 46) |
| 25-29 | 2000 | 116 | 12 (6 to 89) |
| 30-34 | 2000 | 147 | 28 (21 to 108) |
| 35-39 | 2000 | 99 | 61 (48 to 141) |
| 40-44 | 2000 | 83 | 117 (80 to 234) |
| 45-49 | 2000 | 105 | 215 (123 to 318) |
| 50-54 | 2000 | 116 | 296 (175 to 369) |
| 55-59 | 2000 | 112 | 344 (258 to 500) |
| 60-64 | 2000 | 101 | 440 (316 to 457) |
| 65-69 | 2000 | 99 | 462 (321 to 498) |
| 70-89 | 2000 | 371 | 374 (260 to 381) |
| 0-4 | 2001 | 113 | 2 (0 to 23) |
| 5-9 | 2001 | 65 | 1 (0 to 13) |
| 10-14 | 2001 | 51 | 1 (0 to 20) |
| 15-19 | 2001 | 130 | 1 (1 to 46) |
| 20-24 | 2001 | 145 | 5 (2 to 64) |
| 25-29 | 2001 | 127 | 9 (8 to 80) |
| 30-34 | 2001 | 122 | 29 (15 to 112) |
| 35-39 | 2001 | 128 | 45 (37 to 169) |
| 40-44 | 2001 | 107 | 116 (62 to 220) |
| 45-49 | 2001 | 101 | 170 (113 to 306) |
| 50-54 | 2001 | 96 | 281 (165 to 386) |
| 55-59 | 2001 | 91 | 331 (233 to 469) |
| 60-64 | 2001 | 94 | 406 (253 to 482) |
| 65-69 | 2001 | 105 | 472 (298 to 476) |
| 70-89 | 2001 | 391 | 364 (280 to 395) |
| 0-4 | 2002 | 102 | 1 (0 to 16) |
| 5-9 | 2002 | 62 | 0 (0 to 19) |
| 10-14 | 2002 | 64 | 2 (0 to 15) |
| 15-19 | 2002 | 129 | 2 (0 to 45) |
| 20-24 | 2002 | 148 | 8 (3 to 59) |
| 25-29 | 2002 | 120 | 16 (7 to 87) |
| 30-34 | 2002 | 112 | 27 (14 to 116) |
| 35-39 | 2002 | 127 | 62 (33 to 178) |
| 40-44 | 2002 | 98 | 102 (59 to 219) |
| 45-49 | 2002 | 89 | 174 (113 to 307) |
| 50-54 | 2002 | 104 | 241 (176 to 398) |
| 55-59 | 2002 | 116 | 330 (203 to 468) |
| 60-64 | 2002 | 88 | 411 (270 to 452) |
| 65-69 | 2002 | 90 | 435 (315 to 483) |
| 70-89 | 2002 | 384 | 345 (265 to 348) |
| 0-4 | 2003 | 74 | 2 (0 to 20) |
| 5-9 | 2003 | 46 | 1 (1 to 31) |
| 10-14 | 2003 | 59 | 2 (0 to 27) |
| 15-19 | 2003 | 102 | 4 (1 to 44) |
| 20-24 | 2003 | 116 | 6 (4 to 59) |
| 25-29 | 2003 | 137 | 11 (4 to 83) |
| 30-34 | 2003 | 118 | 32 (10 to 123) |
| 35-39 | 2003 | 113 | 53 (27 to 177) |
| 40-44 | 2003 | 109 | 98 (47 to 228) |
| 45-49 | 2003 | 83 | 132 (87 to 304) |
| 50-54 | 2003 | 102 | 230 (126 to 359) |
| 55-59 | 2003 | 92 | 358 (206 to 424) |
| 60-64 | 2003 | 98 | 378 (259 to 456) |
| 65-69 | 2003 | 94 | 418 (297 to 453) |
| 70-89 | 2003 | 342 | 339 (241 to 379) |
| 0-4 | 2004 | 112 | 3 (0 to 28) |
| 5-9 | 2004 | 71 | 2 (0 to 25) |
| 10-14 | 2004 | 81 | 2 (1 to 29) |
| 15-19 | 2004 | 111 | 3 (1 to 54) |
| 20-24 | 2004 | 120 | 9 (2 to 74) |
| 25-29 | 2004 | 136 | 11 (5 to 99) |
| 30-34 | 2004 | 129 | 22 (8 to 126) |
| 35-39 | 2004 | 103 | 50 (23 to 154) |
| 40-44 | 2004 | 128 | 93 (43 to 226) |
| 45-49 | 2004 | 94 | 140 (82 to 240) |
| 50-54 | 2004 | 103 | 208 (135 to 340) |
| 55-59 | 2004 | 98 | 291 (185 to 429) |
| 60-64 | 2004 | 86 | 371 (256 to 414) |
| 65-69 | 2004 | 92 | 414 (311 to 434) |
| 70-89 | 2004 | 312 | 338 (212 to 369) |
plot_state(model, states = c("YearlyAgeInc"), start_time = 59) + theme(legend.position = "none")